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1 .  INTRODUCTION 

The  performance  evaluation  of  finite  element  methods  has  recently  drawn 
a  large  attention.  The  question  is  strictly  related  to  the  development  of  a 
suitable  set  of  benchmark  problems  upon  which  to  verify  and  possibly  to 
compare  the  accuracy  of  different  finite  elements.  In  here  we  refer  to  the 
proposed  set  of  problems  made  by  McNeal  and  Harder  [1],  the  comparison  studies 
of  Robinson  and  Blackham  [2],  [3]  and,  for  a  comprehensive  up-to-date  review 
of  the  matter,  the  proceeding  [4], 

In  this  paper  we  consider  the  benchmark  problem  of  a  simply  supported 
uniformly  loaded  rhombic  plate.  There  are  important  questions  related  to  the 
computation  of  such  a  problem.  For  example,  which  conclusions  can  be 
reasonably  drawn  from  the  results  and  how  to  interpret  the  rapid  deterioration 
of  the  accuracy  as  the  skewness  becomes  larger.  Moreover,  does  the 
computation  of  the  skew  plate  illustrate  adequately  the  sensitivity  of  the 
finite  elements  to  the  skewness  or  is  there  some  other  more  relevant  effect  to 
be  considered?  We  shall  show  that  the  effect  due  to  the  skewness  of  the 
elements  of  the  decomposition  is  negligible,  the  main  effect  being  the 
singularity  of  the  solution  due  to  the  presence  of  obtuse  corners  in  the  plate 
domain. 

The  design  of  benchmark  problems  should  be  made  so  as  to  be  really 
representative  for  a  well  described  class  of  problems  and  effects.  Therefore, 
the  selected  problems  need  to  be  known  in  all  their  aspects  that  can  influence 
the  performance  of  the  finite  element  solution.  For  example,  it  is  well  known 
that  a  special  method  or  approach  can  be  designed  to  perform  very  efficiently 
but  only  for  a  very  narrow  class  of  problems.  Therefore  attention  should  be 
placed  also  on  the  class  of  problems  for  which  the  test  is  characteristic. 
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Otherwise  conclusion  based  on  benchmark  computations  could  be  very 
misleading.  There  are,  of  course,  many  other  questions.  Among  them:  how  to 
design  "academic"  benchmark  problems  isolating  single  effects  which  can  then 
lead  to  useful  conclusions  for  non-academic  environments;  how  to  assess 
robustness  and  reliability  of  the  method;  how  to  estimate  the  computational 
complexity  of  the  method. 

Obviously  conclusions  based  or.  benchmark  computation  can  never  be 
completely  objective.  Nevertheless  useful  information  can  be  drawn  when  based 
on  the  state  of  the  art  of  both  theoretical  and  experience  field.  The 
reliability  has  to  be  understood  not  only  with  respect  to  a  particular 
mathematical  model  but,  in  addition,  also  the  model  has  to  be  considered.  For 
example,  how  accurate  is  the  solution  of  the  Kirchhoff  model  of  a  rhombic 
plate  compared  with  the  solution  of  three  dimensional  linear  elasticity 
problem. 

Throughout  the  paper  we  will  address  the  above  questions,  focusing 
mainly  on  which  type  of  conclusions  can  be  inferred  from  the  results  of 
computation  with  different  finite  elements  of  a  simply  supported  uniformly 
loaded  rhombic  plate.  The  following  aspects  will  be  especially  analyzed: 

1.  Effect  of  the  skewness  of  the  element  of  the  decomposition. 

2.  Effect  of  the  singularity  of  the  exact  solution. 

3.  How  to  improve  the  performance  of  the  finite  element  solution  when 
the  skew  angle  of  the  plate  become  small. 

4.  Relation  between  accuracy  and  computational  complexity  for  various 
finite  element  methods  on  a  given  class  of  problems. 

5.  Class  of  problems  the  benchmark  skew  plate  is  representative  of. 
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The  outline  of  the  paper  is  the  following.  In  Section  2  we  introduce 
the  model  problem  and  characterize,  in  a  suitable  mathematical  way,  the 
properties  of  the  exact  solution  of  the  problem.  Section  3  is  devoted  to  the 
description  of  the  finite  elements  used  for  computations,  together  with  some 
abstract  convergence  results.  In  Section  we  present  the  numerical 
results.  First  we  consider  the  case  of  uniform  decomposition,  then  we- 
consider  an  appropriate  non-uniform  decomposition  allowing  to  considerably 
reduce  the  magnitude  of  the  error.  The  study  of  the  relation  between  accuracy 
and  cost  of  computation  ends  Section  4.  Finally,  in  Section  5  the  conclusions 


are  shown. 
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2.  THE  MODEL  PROBLEM 

We  are  interested  in  the  Kirchhoff  model  for  the  simply  supported  plate 

with  parallelogram  shape.  Let  us  denote  by  ft  the  domain  of  the  (x. ,x?)- 

4 

plane  occupied  by  the  plate,  let  r  =  U  be  its  boundary.  A.,  i  =  1,4, 

i=1 

its  vertices  and  ft  =  ft  U  f.  Let  a  denote  the  skew  angle  of  the  plate. 


Fig.  2.1.  The  domain  of  the  plate. 

We  assume  the  plate  loaded  by  a  transversal  load  g(x1fx2).  We  denote 
by  Hk(Q)  the  standard  Sobolev  space  of  the  functions  with  square-4ntegrable 
derivatives  up  to  tha  order  k.  We  will  also  use  the  Sobolev  spaces  with 
fractional  derivatives  defined  in  the  usual  way  by  interpolation  techniques 
.  (e.g.  the  K-method;  see  [5]  for  details).  We  also  mention  the  Besov  spaces 
Bp  ..(ft)  which  are  very  close  the  the  spaces  Hk(ft),  namely 
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2.1  Hk(G)  c  Bjj>G0(n)  c  Hk'e(Q)  e  >  0 

(see  [5]  for  more  details).  Further  we  denote 

2.2  °H2(ft)  =  {u  €  H2(G):  u|r  -  0}. 

The  exact  solution  Uq  of  our  problem  can  now  be  defined  as  the 

*  ri  p 

minimizer  of  the  quadratical  functional  F(u)  over  H  (Q),  where 


2.3 

2.4  B(u,u) 


2.5 


F( u)  =  B( u, u)  -  Q(u) , 


5  BtJ  t( 


2  2  2 
3  u  +  3U| 

2  2-' 

G  3x^  3x2 


-  2(1-v)( 


2  2 
3  u  3  u 

2  2 
3x^  3x2 


Q(u)  *  {  gCx^Xg)  u  dx-,dx2 

Q 


r  32U 
^  3x^  3x2 


2 

)  )]  dx1dx2} 


2.6 


Et- 


12(1-v2) 


E  is  the  Young's  modulus,  t  the  thickness  of  the  plate  and  v  the 
Poisson's  ratio. 

The  expression  B(u,u)  has  the  meaning  of  the  (strain)  energy  of  the 

plate, 


2.7  Y.|uj  ,  >  |uL  -  B(u,u)y*  >  Y_ ju|  .  ,  Y  >  0, 

:  H  (Q)  n  (Q) 

and  has  all  the  properties  of  a  norm.  Later  we  will  measure  the  error 

2.8  e  *  uq  -  uFE 


between  the  exact  solution  Uq  and  the  finite  element  solution  Upg 


in  the 
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(energy)  norm  defined  by  2.7.  It  is  easy  to  see  that  this  measure  is 
equivalent  to  the  error  in  the  moments  measured  in  L2(fl),  i.e.  in  the  least 
square  way  over  ft. 

Remark  -.1.  The  formulation  2. 3-2. 5  has  proper  meaning  for  a  general  domain 

Q. 

It  follows  easily  by  2.3  and  2.7  that  the  solution  uQ  exists  for  any 
given  load  g(x1 ,x2)  €  H°(Q)  *  L2(fl)  (i.e.  space  of  the  square  integrable 
functions)  and  that  it  is  unique. 

Remark  2.2.  The  solution  exists  for  a  large  class  of  loads.  For  example,  a 
concentrated  load  (g(x1 ,x2)  =  Dirac’s  function)  is  allowed  due  to  the 
inclusion  H2(fi)  C-*.  C°(fl). 

In  the  following  we  will  concentrate  on  the  case  of  loads  given  by 

analytic  functions  on  Q.  A  representative  of  this  class  is  g(x1 ,x2)  =  1 , 

i.e.  uniformly  distributed  load.  (This  example  will  be  considered  in  the 

benchmark  problem).  In  this  case  the  solution  u0  is  an  analytic  function 
4 

on  Q  -  U  Ait  i.e.  uQ  is  not  analytic  at  the  vertices  A^^  but  is 
i=1 

analytic  at  all  the  rest  of  the  boundary.  In  the  neighborhood  of  the  vertex 
A1  (and  A^)  the  solution  has  the  form 

n 

2.9  Un  *  cr*  a  sin  — —  0  +  smoother  terms, 

u  Tr-a 

where  a  (0  <  a  <  ir/2)  is  the  angle  indicated  in  Fig.  2.2a.  In  the 
neighborhood  of  the  vertex  Ajj  (and  A2)  we  have 
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£ 

2.10  un  =  cra  sin  -  e  +  smoother  terms. 

o  a 

With  r,e  we  denote  the  polar  coordinates  with  the  origin  in  A1  (resp.  A2) 
as  shown  in  Fig.  2.2a,b. 


Figure  2.2a,b.  Polar  coordinates  around  the  corners. 

The  solution  uQ  has  a  strongest  singular  behaviour  in.  the  neighbor¬ 
hoods  of  the  vertices  A1  and  A^.  We  see  that  the  singularity  becomes 
stronger  as  a  0  and,  for  any  o,  uQ  f  H^(q).  In  fact,  we  can  show  that, 
for  ~  >  a  >  0 

2.11  u0  €  HT~e(B),  V  e  >  0, 

2.12  f  =  2  +  —  . 

it -a 

Moreover 

2.13  uQ  £  H7(Q)  but  uQ  €  B^,(a). 

Remark  2.3.  The  regularity  of  the  exact  solution  of  the  problem  plays  a 
crucial  role  in  the  performance  of  the  finite  element  method.  In  [6]  we  have 
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characterized  the  smoothness  of  the  solution  in  the  framework  of  countably 
normed  spaces  and  used  those  results  in  relation  with  the  analysis  of  the 
performance  of  the  h-p  version  of  the  finite  element  method  (see  [7,8],  [9], 
[10]). 


Remark  2.4.  The  singularity  of  the  solution  of  the  simply  supported  Kirchhoff 
plate  caused  by  the  corner  of  the  domain  leads  to  some  paradoxical  properties 
of  the  solution.  Consider,  for  example,  the  problem  of  a  simply  supported 
plate  with  Poisson's  ratio  v  =  0,  with  the  shape  of  a  regular  n-sided 
polygon  inscribed  into  a  circle  of  radius  R,  uniformly  loaded  with  a  load 
q.  Let  un  be  the  solution.  Further,  let  u0  be  the  solution  of  the 
analogous  problem  for  the  circular  plate  of  radius  R.  The  solution  un  is 
defined  by  2. 3-2. 5  and  it  is  uniquely  determined  for  every  integer  n,n  =* 
3,4,...  and  un  -*•  u,„  as  n  -*■  ®,  but,  and  this  is  the  paradox,  u,,,  /  uc.  At 
the  center  C  of  the  plate  we  have 


2.-14 


u00(0,0 )  =  lim  un(0,0) 

n->® 


» 


2.15 


uc(0,0) 


f 


where  I  is  the  momentum  of  inertia.  This  means  a  difference  of  more  than 
40J!.  The  implications  of  this  result  (referred  to  as  Babuska  paradox)  have 


been  addressed  in  various  papers  (see  e.g.  [11],  [12],  [13],  [14]). 


3.  THE  FINITE  ELEMENT  APPROXIMATION 

We  are  interested  to  solve  numerically  the  problem  defined  by  2.3~2.5 
using  a  finite  element  approximation.  To  this  end  we  define  the  finite 
dimensional  space  S  of  the  finite  element  solutions  and  then  we  minimize  the 
functional  F(u)  in  2.3  over  the  space  S.  The  core  of  the  finite  element 
method,  in  its  simplest  form,  is  essentially  the  construction  of  a  suitable 
finite  element  space  S.  First  a  triangulation  (or  other  partition)  is 
established  over  the  set  0  and  then  the  space  S  is  constructed.  The 
quality  of  the  finite  element  solution  is  determined  by  the  properties  of  S. 

We  will  consider  in  the  following  three  different  finite  elements  for 
plate  bending  problems: 

a)  Argyris  element  (ARGY) 

b)  reduced  Hsieh-Clough-Tocher  element  (HCTR) 

c)  dual  hybrid  element  (HYBR). 

It  is  well  known  that  a  quite  large  number  of  elements  for  plate  is  described 

in  the  scientific  literature  and  many  of  them  are  implemented  in  the  finite 

element  codes.  A  rough  but  effective  classification  can  be  made  dividing  the 

elements  in  two  categories:  conforming  and  non-conforming,  depending  on 

0  2 

whether  or  not  S  c  VH  (fl).  We  have  restricted  ourselves  to  the  above 
mentioned  elements,  all  of  them  conforming,  although  similar  analysis  as  we 
present  could  be  carried  out  for  other  elements  as  well. 

Let  us  briefly  describe  the  finite  elements  we  have  used,  together  with 


their  basic  properties. 


1  2 


The  Argyris  element 

The  conforrai  ;y  of  a  plate  element  can  be  obtained  in  different  ways.  In 
the  Argyris  element  the  C1 -condition  (i.e.  continuity  of  the  functions  of  S 
and  of  their  first  derivatives),  which  implies  S  c  H  (Q)  is  satisfied 
through  the  use  of  a  complete  space  of  degree  5  (see  [153,  [16]).  This 
corresponds  to  21  degrees  of  freedom  per  triangle,  as  shown  in  Fig.  3.1.  In 
particular  the  value  of  the  displacement,  together  with  its  first  and  second 
derivatives,  is  imposed  at  the  vertices,  the  normal  derivative  is  prescribed 
at  the  midpoint  of  each  side.  We  recall  that  the  number  of  d.o.f.  is  very 
close  of  the  optimal  one  ( 1 8,  see  [173)  needed  to  insure  the  conformity  when 
(only)  a  polynomial  space  is  used. 


Fig.  3.1.  The  Argyris  element  (ARGY). 

The  reduced  Hsieh-Clough-Toeher  element 

The  original  Hsieh-Clough-Tocher  element  is  a  composite  one:  the  basic 
element  is  subdivided  into  three  triangles  and  then  on  each  subtriangle  a 
complete  cubic  polynomial  is  defined.  After  imposing  the  -continuity  at 
the  vertices  and  across  the  sides  of  the  subtriangles  only  12  of  the  initial 
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30  d.o.f.  remain:  displacement  and  its  first  derivatives  at  the  vertices, 
normal  derivative  at  the  midpoint  of  each  side.  The  reduced  element  (see  Pig. 
3.2)  is  obtained  eliminating  the  degrees  of  freedom  normal  derivatives.  This 
correspond  to  require  the  cubic  polynomials  on  each  subtriangle  to  have  linear 
normal  derivatives  on  the  sides.  The  dimension  of  the  space  associated  to  the 
element,  9,  is  optimal  in  the  sense  that  this  is  the  lower  possible  number 
of  d.o.f.  required  for  a  conforming  element.  We  mention  that  this  element  is 
used  in  several  commercial  codes.  For  more  details  we  refer  to  [18],  [19], 
[20]. 


Fig.  3.2.  The  reduced  Hsieh-Clough-Tocher  element  (HCTR). 

The  dual  hybrid  element 

The  basic  idea  for  the  so-called  hybrid  finite  elements  was  first 
suggested  by  Pian  and  Tong  [21].  We  refer  to  [22]  for  a  clear  exposition  of 
the  features  of  this  approach  for  solving  linear  solid  mechanics  problems. 

The  element  we  have  used  has  been  extensively  analyzed,  both  from  theoretical 
'and  numerical  standpoint,  by  Brezzi  and  Marini  (see  [23],  [2h],  [25]).  The 
main  aspect  of  the  finite  element  is  the  particular  choice  of  the  displacement 
approximation  space,  formed  by  cubic  polynomials  only  defined  at  the  boundary 


of  the  triangle.  Obviously,  the  functions  are  in  3ome  way  extended  to  the 
whole  triangle,  but  the  computations  requires  only  the  values  at  the 
interelement  boundaries  The  degrees  of  freedom  are  shown  in  Fig.  3*3  and 
consists  of  values  of  displacement  and  its  first  derivatives  at  the 
vertices.  The  dimension  of  the  approximation  space  is  9. 


Fig.  3.3.  The  dual  hybrid  element  (HYBR). 


Let  us  now  recall  the  approximation  properties  of  the  elements  we  have: 

THEOREM  3.1.  Let  the  triangulation  of  8  be  quasi-uniform (satisfying  the 
minimal  angle  condition),  the  solution  Ug  €  H^r)  fl  ®H^(8),  k  >  2  (integer 
or  fractional).  Then the  following  estimate  holds: 


3.1  lU0  UFE^E  -  Ch^U°^Hk(Q) 
with 

3.2  u  a  min(k-2,4)  for  ARGY 

3.3  \i  =  min(k-2,1)  for  HCTR,  HYBR, 


where  h  is  the  maximum  size  of  the  elements  of  the  decomposition,  C  is  a 
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constant  depending  on  the  aspect  ratio,  the  skewness  of  the  elements  and  the 
type  of  the  finite  element. 

Theorem  3.1  is  related  to  the  standard  h-version  of  the  finite  element 
method.  As  the  number  K  of  degrees  of  freedom  is  (asymptotically)  of 
order  we  can  re-write  relation  3.1  in  the  form 

l-*  lu0  '  uFe!e  S  k,  • 

Theorem  3*1  is  formulated  in  a  general  way.  For  the  problem  of  uniformly 
loaded  simply  supported  skew  plate  we  can  state  the  following: 

THEOREM  3,2.  Let  the  triangulation  of  fl  be  quasi -uniform  (satisfying  the 
minimal  angle  condition)  and  the  load  be  .uniform  on  R.  Then  there  exist  two 
constant  and  C2  such  that  the  following  estimate  holds: 

-*/a  -2-  -4-2- 

3.5  C,H  *  “  <  Iu0  -  upE|E  <  0,H  ’  “  , 

where  the  constants  C1  and  C-2  depend  on  the  aspect  ration,  the  skewness  of 
the  elements,  the  finite  element  itself  (ARGY,  HCTR  or  HYBR) ,  but  are 
independent  of  N  (i.e.  the  rate  of  convergence  is  the  same  for  all  three 
elements ) . 

Theorem  3-2  shows  that  the  rate  of  convergence,  and  therefore  the 
accuracy,  deteriorates  while  the  skew  angle  o  decreases.  This  behaviour  is 
almost  independent  of  aspect  ratio  and  skewness  of  the  elements  as  these 
influence  only  the  constants  and  C2.  The  theorem  holds  not  only  for  our 
choice  of  elements  but  it  is  more  general,  the  only  exception  being  the  case 
when  special  elements  with  the  shape  functions  of  the  form  2.9  are  used  (this 
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is,  for  example,  the  reason  of  the  performance  of  the  element  ELFIN  in  the 
comparative  study  made  by  Robinson  [26]). 

Theorem  3.2  has  only  asymptotic  character  when  the  constants  are  not 
specified.  Obviously,  in  practice,  these  constants  are  important  and  benchmark 
computations  may  characterize  them  well  for  particular  classes  of  problems. 

£o  far  we  ha va  addressed  only  the  case  of  quasi-uniform  mesh.  If  the 
mesh  is  properly  refined  then  we  have  the  following: 

THEOREM  3.3.  Let  the  load  acting  on  the  plate  be  uniform  on  Q  and  the  mesh 
decomposition  properly  selected.  Then  there  exist  constant  such  that  the 

following  estimate  holds: 

3.6  |u0  -  uFElE  <  C,H-^  . 

with  n  =  4  for  ARGY^  \\  -  1  for  HCTR  and  HYBR.  Further  for  any  mesh 
C2N"V^  <  |  Uq  -  uFE|E  where  C2  >  0  depends  only  on  3kewnees  of  the  elements.. 

D 

Theorem  3.3  shows  that  the  performance  of  ARGY  is  especially  good  when 
proper  design  of  the  mesh  is  made.  We  will  show  that  this  choice  will  permit 
the  acnievement  of  a  higher  accuracy. 

Theorems  3.1,  3*2  and  3-3  follow  from  the  standard  mathematical  theory 
of  finite  elements  (see  e.g.  [19],  [27]). 

So  far  we  have  discussed  only  the  performance  with  respect  to  the  energy 
norm  measure  of  the  error.  Analogous  behaviour  is  expected  for  other  norms  of 
the  accuracy.  We  will  partially  address  questions  of  this  type  in  the 


following  section 
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4.  THE  BENCHMARK  PROBLEM  COMPUTATION 

In  this  section  we  compare  the  performance  of  the  finite  elements 
described  in  Section  3.  We  will  analyze  energy  and  displacement  errors  with 
respect  to  the  numbers  of  degrees  of  freedom  and  the  cost  of  computational 
work.  Quite  often  in  the  literature  only  the  relation  between  error  and 
degrees  of  freedom  is  shown.  Although  this  is  a  very  important  characteriza¬ 
tion,  the  most  important  information  is  the  relation  between  accuracy  and 
computational  cost.  The  latter  depends,  of  course,  on  different  factors 
(programming  technique,  solver  algorithm,  etc.)  which  make  it  not  completely 
well  defined  while  the  first  characterization  has  a  completely  precise 
meaning.  Nevertheless,  using  standard  level  of  programming  and  suitable 
description  of  the  machine  cost  very  reliable  informations  can  be  obtained. 


4.1.  THE  UNIFORM  MESH  DECOMPOSITION 

We  first  consider  the  case  of  uniform  mesh.  In  Fig.  4.1  an  example  of 
such  a  mesh  is  shown. 

We  have  solved  the  skew  plate  problem  for  four  different  values  of  the 
skew  angle  a:  80°,  60°,  40°,  30°. 


Fig.  4.1.  A  uniform  mesh  with  triangular  elements. 
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The  physical  data  of  the  plate  are  the  following: 


side  length 

1.0 

in 

Young's  modulus 

3.  *  107 

lb/ in2 

Poisson's  ratio 

0.3 

thickness 

0.1 

in 

load 

1. 

lb/in2 

First  we  show  the  relation  between  the  energy  norm  of  the  error  and  the 
number  of  degrees  of  freedom.  Let  EEX  denote  the  exact  energy  of  the 
plate,  EpE  the  energy  of  the  finite  element  solution.  The  relative  energy 
norm  |eIER  fche  error  e  =  uQ  -  uFE  can  be  expressed  (using  basic 
properties  of  the  finite  element  solution)  in  the  following  way 
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Figures  4.2a-d  show  the  results  for  different  values  of  a  when  a  mesh  of  the 
type  shown  in  Fig.  4.1  is  used.  In  each  figure  results  for  ARGY,  HCTR  and 
HYBR  elements  are  given.  Both  the  scales  are  of  logarithmic  type.  In  all  the 
cases  the  ARGY  element  give  better  performances.  This  is  expected  when  the 
singularity  is  still  weak  (i.e.  a  *  80°),  but  it  happens  for  each  value  of 
a  which  supplies  the  information  about  the  constants  C1  and  C2  (see 
Theorem  3.2).  Within  each  test  the  order  of  convergence  is  nearly  the  same 
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Fig.*1.2a-d.  Energy  norm  of  the  error  against  number  of  degrees  of  freedom  for 
uniform  mesh.  The  slope  of  the  triangle  denotes,  in  each  figure,  the  theoret¬ 
ical  rate  of  convergence. 
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for  the  three  different  methods.  In  the  figures  the  theoretical  order  is 
shown.  As  the  error  is  (see  3.5}  of  order 

-y  _2_ 

4.2  N  1T“° 

the  slope  of  the  lines  is  expected  to  be  close  to  the  value  -V2  — ^  ♦  A  quite 
good  agreement  between  theoretical  and  computational  order  of  convergence  can 
be  seen.  We  note,  as  expecteu,  the  rapid  increasing  of  the  error  when  a 
becomes  smaller.  In  particular  in  Fig.  4.2a,b,c  the  same  scales  have  been 
used  to  emphasize  the  deterioration  of  the  accuracy  while  a  assumes  the 
values  80,  60,  40.  Due  to  the  magnitude  of  the  error  a  different  scale  is 
used  in  Fig.  4. 2d  (a  «  30°). 

Let  us  now  consider  the  displacement  of  the  plate.  Let  uM0(°)  denote 
the  value  at  the  center  of  the  plate  computed  by  Morley  ([28],  [29])  using 
series  expansion,  uFE(C)  the  finite  element  solution.  The  relative 
displacement  error  is  simply  defined  as 


4.3 


M  r  M0(C)  UFE(C\ 

D*  =  l - - (rfS - J 


100. 


Figures  4.3a-d  give  the  error  D%  against  number  of  degrees  of  freedom.  The 
behaviour  of  the  displacement  error  is  roughly  the  same  as  the  error  in  the 
energy^  aeeau  This  is  not  surprising  due  to  the  relation  linking  energy  and 
displacement  (i.e.  the  error  in  the  energy  is  the  average  error  in 
displacement). 

Let  us  now  explain  the  reason  of  the  deterioration  of  accuracy  as  a 
becomes  smaller.  Obviously  when  a  decreases  the  skewness  of  the  elements 
increases  and,  at  the  same  time,  the  results  worsen.  This  fact  sometimes 
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Fig.  M.3a-d.  Relative  displacement  error  at  the  center  of  the  plate  against 
number  of  degrees  of  freedom  for  uniform  mesh.  The  slope  of  the  triangle 
denotes,  in  each  figure,  the  theoretical  rate  of  convergence. 
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seems  to  lead  to  the  conclusion  (see  [26])  that  the  effect  of  the  skewness  of 
the  elements  is  the  reason  of  worsening.  The  major  factor  is  the  change  of 
the  smoothness  of  the  exact  solution  in  dependance  of  a  and  not  the  skewness 
of  the  elements.  To  illustrate  this  fact  let  us  compute  the  error  when  meshes 
of  type  shown  in  Fig.  4.4  b  and  d  (in  Fig.  4.4  a  and  c  the  usual  meshes  are 
shown).  Of  course.,  the  skewness  of  meshes  4.4  b  and  d  is  larger  than  the  one 
of  4.4  a  and  c,  respectively.  Figures  4.5a-d  compare  the  measure  in  the 
energy  norm  of  the  error  for  both  normal  and  'skewed'  meshes.  The  continuous 
lines  refer  to  the  normal  mesh,  the  dashed  lines  to  the  'skewed'  mesh. 

Figures  4.6a~d  show  analogous  results  for  the  displacement  error. 


-Fig.  4.4a-d.  Examples  of  mesh  used  to  evaluate  the  sensitiveness  of  the 
numerical  solution  to  the  skewness  of  the  elements  of  the  decomposition. 
Left:  'natural'  mesh;  right:  'skewed'  mesh. 
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We  see  that  the  singularity  of  the  solution  and  not  the  skewness  of  the 
elements  is  the  main  reason  for  the  change  of  performance  of  the  method.  Only 
a  relatively  slight  contribution  is  given  by  the  'skewed'  mesh  negligible 
compared  with  the  influence  of  the  singularity. 

Remark  4.1.  Of  course  the  form  of  the  triangle,  more  precisely  its  maximal 
angle  (see  [30])  has  essential  influence  on  convergence  (i.e.  when  it  is 
violated  convergence  does  not  occur  at  all).  Nevertheless,  in  our*  case  the 
maximal  angle  of  any  triangle  is  far  enough  from  ir  and  thus  its  effect  is 
negligible  compared  with  the  smoothness  effect. 

4.2.  The  non-uniform  mesh  decomposition 

The  presence  of  singular  behaviour  at  the  obtuse  corners  of  the  plate 
makes  the  moments  unbounded.  The  rate  of .convergence  (see  4.2)  is  so  small 
that  uniform  mesh  gives  completely  unacceptable  results.  In  example,  for 
a  *  30°,  extrapolating  the  error  behaviour  shown  in  Fig.  4.3d,  we  see  that  a 
number  of  degrees  of  freedom  N  >  10®(1!)  is  needed  to  achieve  a  5%  accuracy 
for  the  displacement  at  the  center  of  the  plate.  If  we  properly  refine  mesh 
then  Theorem  3*3  applies  and  much  better  results  can  be  achieved,  especially 
for  the  ARGY  element.  The  optimal  meshes  need  to  be  graded  in  the  place  where 
singular  solution  occurs.  Figures  4.7a-d  show  a  sequences  of  refined  meshes 
used  in  the  computation.  Figures  4,8a-d  show  the  error  in  the  energy  norm 
against  the  number  of  degrees  of  freedom  for  angles  a  =  80°,  60°,  40°,  30°. 
The  performance  of  the  ARGY  element  with  the  non-uniform  meshes  of  Fig.  4.7  is 
reported,  together  with  the  results  related  to  uniform  meshes.  Also  the  slope 
related  to  the  theoretical  order  given  in  Theorem  3.3  is  shown  in  the  figures. 
The  improvement  of  the  accuracy  Is  really  effective,  especially  when  o 


becomes  smaller 
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The  mesh  shown  in  Fig.  4.7  are  geometric  meshes  with  a  factor  V2.  This 
type  of  decomposition  is  effective  as  far  as  the  major  contribution  to  the 
error  appears  in  the  smallest  element,  due  to  the  singularity  of  the  solution. 
If  the  singularity  is  not  strong  enough  such  a  geometric  mesh  is  not  effective 
because  the  major  error  comes  out  of  the  larger  internal  element.  This 
behavior  can  be  seen  in  Fig.  4.8a  (a  =  80°,  weak  singularity).  However, 
after  a  certain  number  of  successive  refinements  analogous  situation  will 
occur  for  the  other  angles  as  well.  When  this  phase  is  reached  simultaneous 
refinement  around  the  singularity  location  and  outside  should  be  made. 

Figures  4.9a-d  show  the  error  in  the  displacement  at  the  center  of  the 
plate.  Again,  the  slope  of  the  theoretical  convergence  rate  (that  is  twice 
the  order  of  the  energy  norm)  is  indicated. 

We  see  that  the  ARGY  element,  which  is  of  higher  degree,  performs  very 
well  in  the  case  of  non-uniform  mesh.  This  is  among  others  due  to  the  fact 
that  higher  order  finite  elements  are  able  to  better  "absorb"  the  singulari¬ 
ties  of  the  solution,  as  shown  in  [32],  [33]  in  relation  to  the  p-version  of 
the  finite  element  method. 

It  is  necessary  to  note  that  ARGY  element  is  over cons trained  (i.e.  the 
second  derivatives  at  the  vertices  are  used).  This  can  sometimes  lead  to 
difficulties  (not  present  in  the  problem  we  are  dealing  with).  For  example, 
when  the  plate  change  thickness  the  second  derivatives  of  the  solution  a^e  not 
continuous  while  the  overconstraining  enforces  continuity.  Nevertheless,  a 
simple  adaptation  of  the  elements,  using  the  fact  that  the  character  of  the 
discontinuity  is  known,  can  avoid  the  difficulty.  Another  well  known  cases 
are  those  when  boundary  condition  changes  from  clamped  type  to  nonclamped  and 
in  corners  of  clamped  boundary  (which  can  be  delt  with  by  a  proper  refinement 
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Fig.  4.9a-d.  Relative  displacement  error  at  the  center  of  the  plate  against 
number  of  degrees  of  freedom  for  uniform  and  non-uniform  mesh. 
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of  the  mesh).  Hence  the  results  of  the  discussed  benchmark  is  representation 
for  simply  supported  polygonal  plates. 

4.3.  Accuracy  and  computational  work 

The  comparison  accuracy  -  number  of  degrees  of  freedom  is  the  most  used 
in  the  scientific  literature.  However,  given  a  problem  to  solve  within  a 
certain  range  of  precision,  the  most  important  relation  is  between  accuracy 
and  cost  of  computation.  This  should  be,  and  indeed  it  is  in  practical 
environment,  the  criterium  upon  which  to  rest  the  selection  of  the  most  suit¬ 
able  finite  element.  Only  recently  this  problem  has  been  approached  from  a 
quite  general  standpoint  (see  [34]).  Hereafter  we  will  give  only  some  general 
lines  of  this  important  question,  focusing  mainly  on  the  numerical  results 
obtained  with  the  finite  elements  previously  described.  The  finite  element 
computation  requires  the  execution  of  the  following  steps: 
a  -  topology,  mesh  generation 
b  -  local  stiffness  matrices 
c  -  assembly  and  solution 

d  -  postprocessing  (i.e.  computation  of  required  data). 

Sometimes  assembly  is  combined  with  computation  of  local  stiffness  matrices. 
The  cost  of  the  mesh  generation  and  the  cost  of  postprocessing  are  not  too 
much  sensitive  to  the  elements.  On  the  other  hand,  parts  b  and  c  depend 
heavily  on  the  type  of  finite  element.  Toward  a  simple  approach  we  will  only 
consider  the  cost  of  parts  b  and  c.  As  regards  the  accuracy,  we  have  already 
seen  that  different  measures  can  be  used  (e.g.  energy  norm,  L„  norm,  etc.). 

We  will  especially  consider  the  energy  norm  although  pointwise  accuracy  could 
be  examined  as  well.  The  conclusion  we  will  base  on  the  energy  norm  will 
essentially  hold  for  pointwise  accuracy. 


Computer  time  depends  on  both  hardware  (operative  system,  compiler, 
etc.)  and  software  (programmation  technique,  etc.).  Nevertheless,  it  is 
reasonable  to  assume  that  the  state  of  the  art  programming  will  lead  to 
acceptable  and  objective  conclusions.  Of  course,  some  very  special  techniques 
could  influence  the  conclusions  but  this  has  to  be  considered  as  a  "different" 
method  (as  it  is  the  case  for  parallel  computer  oriented  methods,  which  are 
not  optimal  for  sequential  machines  and  could  give  better  performances). 

The  results  presented  throughout  the  paper  has  been  obtained  using  an 
Apollo  410  computer.  The  programs  are  performed  in  a  standard  way  and  are 
included  in  the  finite  element  library  MODULEF  (see  [35]).  The  final  linear 
system  is  solved  via  elimination  method  based  on  the  Choleski  factorization. 

In  Fig.  4.10a,b  the  time  for  the  computation  of  local  stiffness 
matrices,  resp.  for  assembly  and  solution,  against  the  number  of  degrees  of 
freedom  is  shown.  As  expected,  in  both  the  figures  the  time  spent  is  larger 
for  ARGY  than  HCTR  and  HYBR  element,  depending  mainly  on  the  dimension  of  the 
elementary  stiffness  matrix. 

As  previously  said  we  are  interested  in  accuracy  versus  cost  and  this  is 
the  subject  of  Fig.  4.11a-d  where  the  cost  is  represented  as  total  time  for 
computation  of  local  stiffness  matrices,  assembly  and  solution  while  the  error 
in  tne  energy  norm  is  chosen  for  the  accuracy. 

Analogous  result  hold  for  the  displacement  error  at  the  center  of  the 
plate.  We  compared  the  performance  of  the  elements  ARGY,  HCTR  and  HYBR  on  a 
uniform  mesh.  In  addition,  we  have  also  shown  the  performance  of  the  ARGY 
element  for  refined  mesh.  (Such  a  mesh  needs  more  work  to  generate.)  The 
proper  mesh  refinement  has  larger  effects  on  the  performance  of  the  higher 
degree  element  than  the  lower  degree,  hence,  for  the  refined  mesh  the 
comparison  will  be- still  more  favorable  for  ARGY. 
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Fig.  4.10a,b.  Time  for  computation  of  local  stiffness  matrices  against  number 
of  degrees  of  freedom  (a).  Time  for  assembly  and  solution  against  number  of 
degrees  of  freedom  (B)  (HCTR  and  HYBR  connect). 

Remark  4.1 .  We  analyzed  the  performance  of  the  finite  element  method  to  solve 
the  plate  problem  according  to  the  mathematical  model  based  upon  the  Kirchhoff 
theory.  Of  course,  the  question  about  the  validity  of  such  a  model  is  essen¬ 
tial.  We  can  formulate  the  plate  problem  as  a  3-dimensional  elasticity  pro¬ 
blem  and  give  it  (exact)  solution  the  meaning  of  true  solution.  Then  we  can 
compare  this  value  with  the  (exact)  solution  of  the  Kirchhoff  (2-dimensional) 
model. 

It  is  well  known  that  as  the  thickness  t  0  the  3-dimensional 
solution  converges  to  the  2-diraensional  one  (see,  e.g.  [36],  [37],  [38], 

CM],  [41]).  The  accuracy  of  the  Kirchhoff  solution  depends  very  much  on  both 
thickness  t  and  skew  angle  a. 
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Denote  by  u,  v,  w  the  3D  displacement  of  the  plate  and  let  us  model  the 
simple  support  by  constraining  w  -  0  at  the  sides  boundary  surface  while 
leaving  u  and  v  unconstrained.  Let  E^q  be  the  exact  (strain)  energy  of 
the  3D  formulation,  EK  the  Kirchhoff  one.  Then  the  quantity 

"  Ev  ‘4 

H.4  EE  =  (-2| - -) 

E3D 

is  very  close  to  the  relative  error  of  the  Kirchhoff  solution  in  the  energy 
norm  (like  in  relation  4.1).  In  fact,  for  v  ®  0,  is  the  exact  relative 

=  1 ,  we  have : 

j  WoD  dx j  dx2, 

Q 

J  wK  dx,  dx2, 
a 

J  (w3D~wK)dx1  dx2» 
a 

where  wgD,  wk  denote  the  vertical  displacement  in  the  3D  and  Kirchhoff 
model,  reap.  Hence,  ED  *  (EE)^  is  the  relative  average  error  in  the  dis¬ 
placement  at  the  loaded  surface. 

Table  4.1  shows  the  value  EE,  ED  for  t  =  .0.01  and  various  values 
a.  The  values  EgD  have  been  computed  out  of  the  results  obtained  using  the 
code  STRIPE  (with  p-version  capability),  developed  by  the  Computational 
Mechanics  Center  of  the  Aeronautical  Research  Institute  of  Sweden  (see  [39]), 
on  the  CRAY  X-MP/48  of  the  Pittsburgh  Supercomputing  Center  (PSC).  A  careful 
extrapolation  technique,  based  on  the  results  for  different  values  of  p  and 
meshes,  has  been  applied.  The  values  EK  have  been  computed  with  analogous 


error  measured  in  the  energy  norm. 
Since  the  uniform  load  is  g 

4.5  E3D  = 

4.6  Ek  - 

4.7  ®3D  ”  EK  = 


extrapolation  technique,  based  on  the  results  for  different  values  of  p  and 
meshes,  has  been  applied.  The  values  EK  have  been  computed  with  analogous 
procedures  out  of  the  numerical  results  described  in  the  subsections  4. 1-4. 3. 

TABLE  4.1 
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30 

38.17 

14.57 

ERR1  :  %  ERROR 

IN  THE 

ENERGY  NORM 

ERR2 :  t  AVERAGE  DISPLACEMENT  ERROR 

thickness  -  0.01 
Poisson's  ratio  =0.3 

We  see  that  the  Kirchhoff  model  gives  unacceptable  results  except  for 
plates  either  close  to  the  square  shape  or  with  very  small  thickness.  In  a 
forthcoming  paper  we  will  address  in  details  the  reliability  of  various  plate 
models  of  Reissner-Mindlin  type  and  discuss  the  accuracy  of  their  finite 
element  approximation. 
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5.  CONCLUSIONS 


Let  us  summarize  the  conclusions  which  can  be  drawn  from  the  benchmark 
computation  of  the  uniformly  loaded  simply  supported  skew  plate  with  Kirchhoff 
formulation. 

1.  The  major  reason  of  the  error  dependance  on  the  skewness  of  the 
plate  is  the  singularity  of  the  solution  and  not  the  skewness  of  the  element 
of  the  decomposition.  Hence  this  benchmark  problem  is  not  well  suited  for 
comparing  different  elements  except  for  the  case  when  the  performance  in 
presence  of  strong  singularity  of  the  solution  needs  to  be  evaluated. 

2.  The  problem  characterizes  well  any  polygonal  simply  supported  plate 

subject  to  analytic  load.  The  solution  in  the  neighborhood  of  a  critical 

vertex  with  internal  angle  3  (see  Fig.  5.1)  has,  for  any  3  such  that 
-1 

(«■)  i3  not  integer,  the  following  expression: 

irK.  irKp 

2+  *Ki  ~~  ttK 

u  =  c-]r  sin  -j-  0  +  C2r  p  sin  e  +  smoother  terms 


with  K-, , 
ttK? 

and  —  > 


1 . 


integers  and  since  u  €  H2(a)  we  have  to  have 
Hence  the  major  singularity  term  is: 


ttK, 


2  + 


>  1 


u  =  CrS  sin  j  6  for  0  <  3  <  it 

2"  B  TT 

u  -  Cr  B  sin  ~  0  for  it  <  3  <  2tt. 

For  example,  the  L~shaped  domain  with  3  =  ^  tt  should  be  compared  with  a  skew 
plate  with  skewness  angle  a  =  45°. 
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Fig.  5.1.  A  reentrant  corner  in  a  L-shaped  domain. 

Therefore,  with  proper  correspondence  the  benchmark  problem  can  be  applied 
quite  generally. 

3.  The  higher  order  elements  provide  better  accuracy  (in  the  engineer¬ 
ing  range)  than  the  lower  ones  for  the  same  computational  cost,  both  in  pres¬ 
ence  of  singularity  and  for  smooth  solution,  either  for  uniform  or  properly 
designed  meshes.  By  accuracy  we  mean  the  error  in  the  energy  norm  or  point- 
wise  error.  The  computer  cost  takes  into  account  computation  of  local  stiff¬ 
ness  matrices,  assembling  procedure  and  (direct)  solving  time.  The  relation 
between  accuracy  and  cost  are  in  agreement  with  the  theoretical  model 
introduced  in  [3*0. 

**.  Because  anisotropic  plates,  after  proper  coordinates  transformation, 
has  the  same  (or  very  similar)  properties  of  isotropic  ones,  the  conclusions 
-hold  also  for  this  case.  Of  course,  the  notion  of  the  angle  is  now  not  only 
geometrical  but  depends  also  on  the  anisotropy.  For  example,  the  finite  ele- 
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ment  of  an  anisotropic  square  plate  behaves  then  as  the  solution  of  an 
isotropic  skew  plate. 

5.  Accuracy  of  the  Kirchhoff  model  compared  with  3  dimensional  linear 
elasticity  solution  could  not  be  acceptable  for  larger  skewness  of  the  plate. 
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